3. DADES DE MENORCA
Maria Àngels Jaume
2023-03-07
1. LECTURA DE DADES
En aquest document farem el mateix estudi que amb les dades de Mallorca però utilitzant les de Menorca. Els càlculs i procediments són anàlegs al de l’illa de Mallorca, per això obviarem alguns detalls.
En primer lloc, grafiquem en forma de sèrie temporal aquestes dades:
turisme <- read_excel("IBESTAT.xls")
turisme$Data <- gsub("M","-",turisme$Data)
gastos.ts<-ts(turisme[-1], start=c(2015,10), frequency = 12)
men <- data.frame(x=1:86, y=turisme$Men) #dades de Menorca
serie_men <- ts(men$y,start=c(2015,10),frequency = 12)
plot.ts(serie_men, main="Menorca", xlab="Any", ylab="Despeses mensuals en €")
Les dades van des d’octubre de 2015 a novembre de 2022 (llavors tenim 7 cicles complets). Podem observar que l’efecte de la COVID és veu clarament en 2020.
Per veure millor l’estacionalitat visualitzem cada període mensualment:
seasonplot(serie_men, col=c("brown", "blue","red", "orange", "pink", "purple", "yellow","green"),year.labels=TRUE, main="Estacionalitat de Menorca", xlab="Mes", ylab=" Despeses en €")
legend("bottomright",
legend = c(2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022),
fill = c("brown", "blue","red", "orange", "pink", "purple", "yellow","green"),
border = "black")
Observam que efectivament, així com es comenta en el treball, el període de confinament va ser de principi de març, on comença a decréixer el valor dels preus, fins a finals de juny, on pareix que s’han tornat a recuperar els valors anteriors a la COVID.
Dividirem la sèrie en tres trams:
Des de 10-2015 fins 12-2018 (serie1_men), amb el que predirem un període. És a dir predirem de 1-2019 a 12-2019. (amb aquest comprovam que els mètodes funcionen)
Des de 10-2015 fins 12-2019 (serie2_men), amb la que predirem el període de la COVID (és a dir, l’any 2020)
Des de 10-2015 fins 12-2020 (serie3_men), amb la que predirem el períodes després de la COVID (és a dir, l’any 2021)
Visualitzem-los:
serie1_men <- ts(men$y[1:39],start=c(2015,10),frequency = 12)
serie2_men <- ts(men$y[1:51],start=c(2015,10),frequency = 12)
serie3_men <- ts(men$y[1:63],start=c(2015,10),frequency = 12)
par(mfrow=c(2,2), mar=c(4,4,4,1)+.1)
plot.ts(serie1_men, main="Serie abans de la COVID \nmenys un cicle", xlab="Any", ylab="Euros")
plot.ts(serie2_men, main="Serie abans de la COVID", xlab="Any", ylab="Euros")
plot.ts(serie3_men, main="Serie abans i durant de la COVID", xlab="Any", ylab="Euros")
plot.ts(serie_men, main="Serie completa", xlab="Any", ylab="Euros")
2. AJUST DELS MODELS:
Abans d’aplicar algun model, estudiem un poc la sèrie:
serie_men.lm <- lm(y~x, data=men)
plot.ts(men$y, main = "Menorca", ylab = "Despeses mensuals en €", xlab="Índex de cada mes")
#dibuixam la recta de regressió sobre els punts
abline(serie_men.lm, col='red')
Si dibuixam la recta de regressió sobre les nostres dades, tot i que aquestes estan molt disperses i no s’ajusten bé a la recta, podem observar que la tendència és una mica creixent, encara que és manté més o menys constant.
boxplot(serie_men~cycle(serie_men), xlab = "mesos", ylab = "Despeses en €", main="Boxplot de Menorca")
Podem observar també la presència d’estacionalitat, que pren els seus valors màxims a la temporada d’estiu i els seus mínims en l’hivern, fet que corresponen amb les dades turístiques a les illes.
Aplicarem diversos models per ajustar la nostra sèrie i fer-ne una predicció per llavors comparar quin és el millor.
Vegem com actua cada model:
2.1. ABANS DE LA COVID
- REGRESSIÓ SEGMENTADA
Com hem comentat anteriorment la recta de regressió no s’ajusta bé a les dades, de fet el valor del \(R^2\) és molt baix:
serie1_men.lm <- lm(y~x, data=men[1:39,])
summary(serie1_men.lm)$adj.r.squared## [1] 0.03739615Per això utilitzam el paquet segmented per ajustar a una
recta de regressió a trossos.
Anem a utilitzar la funció selgmented() per veure quants
de punts de canvi detecta:
punts_canvi_serie1_men <-selgmented(serie1_men.lm, Kmax = 10, type="bic")## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 .. 7 .. 8 .. 9 .. 10 ..
##
## BIC to detect no. of breakpoints:
## 0 1 2 3 4 5 6 7
## 423.9598 430.9145 432.3182 431.8441 426.2749 397.3849 386.0818 390.6075
## 8 9 <NA>
## 392.2832 399.0984 400.0984
##
## No. of selected breakpoints: 6Aplicam la funció segmented() que ens calcula la
regressió segmentada:
serie1_men.seg <- segmented(serie1_men.lm, seg.Z=~x, psi = c(4,12,16,23,28,35))Aquesta funció ens demana que introduïm els valors on es troben els punts de canvi, i aquesta ens dona el valor estimat. Els valors reals dels punts de canvi son:
summary(serie1_men.seg)$psi## Initial Est. St.Err
## psi1.x 4 4.206349 0.6612491
## psi2.x 12 11.801361 0.3959586
## psi3.x 16 15.390922 0.3696944
## psi4.x 23 22.697828 0.4490507
## psi5.x 28 28.566720 0.4342304
## psi6.x 35 35.595164 0.3499531Que corresponen, aproximadament, a: 1-2016, 9-2016, 12-2016, 8-2017, 2-2018 i 9-2018.
Obtenim que el valor de \(R^2\) és bastant alt
summary(serie1_men.seg)$adj.r.squared## [1] 0.8366384Anem a visualitzar la regressió segmentada damunt les nostres dades
p_serie1_men <- ggplot(men[1:39,], aes(x=x, y=y)) + geom_line()+
ylim(c(0,1300))+
ggtitle('Regressió lineal i segmentada sobre les dades \nde Menorca abans de la COVID') +
xlab('Índex del mes') +
ylab('Despeses mensuals en €')
my.coef1_men <- coef(serie1_men.lm) #coeficients de la recta de regressió lineal
p_serie1_men <- p_serie1_men + geom_abline(intercept=my.coef1_men[1], slope=my.coef1_men[2], color="green") #afegim la recta de regressió lineal
my.model1_men <- data.frame(x=1:39, y=fitted(serie1_men.seg)) #model nou amb els valors ajustats de la regressió segmentada
p_serie1_men <- p_serie1_men + geom_line(data=my.model1_men, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada
my.lines1_men <- serie1_men.seg$psi[,2]# guardam on estan els punts de canvi estimats
p_serie1_men <- p_serie1_men + geom_vline(xintercept=my.lines1_men, linetype="dashed")
p_serie1_men <- p_serie1_men + theme(panel.background = element_blank()) + #eliminam el fons i la quadrícula del gràfic
geom_vline(xintercept=0) + geom_hline(yintercept=0)
ggplotly(p_serie1_men)Visualment també es pot observar que la recta de regressió segmentada s’ajusta millor a les nostres dades.
Anem a calcular les equacions d’aquestes rectes, sabem que les rectes tenen la forma \(y=mx+n\), on \(m\) és la pendent i \(n\) el valor de tall de l’eix y.
Hi ha una funció del paquet segmented que ens dona
aquesta valors de \(n\):
intercept(serie1_men.seg)## $x
## Est.
## intercept1 824.68
## intercept2 139.10
## intercept3 3347.90
## intercept4 -1263.70
## intercept5 3720.40
## intercept6 -2821.90
## intercept7 8470.20També tenim una altra funció que ens calcula les pendents:
slope(serie1_men.seg)## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -87.227 39.496 -2.2085 -168.570 -5.8827
## slope2 75.759 16.690 4.5391 41.385 110.1300
## slope3 -196.140 39.496 -4.9660 -277.480 -114.7900
## slope4 103.490 16.690 6.2008 69.119 137.8700
## slope5 -116.090 21.112 -5.4990 -159.570 -72.6130
## slope6 112.920 16.690 6.7658 78.549 147.3000
## slope7 -204.310 39.496 -5.1729 -285.660 -122.9700Aleshores, la nostra recta segmentada és:
\(\hat{y}= \left\{ \begin{array}{lcc} -87.227x + 824.68 & si & x \leq \psi_1 \\ \\ 75.759x + 139.10 & si & \psi_1 < x \leq \psi_2 \\ \\ -196.140x + 3347.90 & si & \psi_2 < x \leq \psi_3 \\ \\ 103.490x - 1263.70 & si & \psi_3 < x \leq \psi_4 \\ \\ -116.090x + 3720.40 & si & \psi_4 < x \leq \psi_5 \\ \\ 112.920x - 2821.90 & si & \psi_5 < x \leq \psi_6 \\ \\ -204.310x + 8470.20 & si & \psi_6 < x \\ \end{array} \right.\)
on \(\psi_1= 4.21\), \(\psi_2=11.8\), $_3=15.39 $, $_4=22.7 $, $_5=28.57 $ i \(\psi_6= 35.6\)
Vegem els errors:
accuracy(serie1_men.seg)## ME RMSE MAE MPE MAPE MASE
## Training set -2.915759e-15 70.70987 57.10189 -1.192331 8.628459 0.2902722Anem a fer la predicció de 1-2019 a 12-2019.
El darrer punt de canvi que tenim és en setembre de 2018 i sabem, seguint el patró obtingut amb les dades d’entrenament, que els següents és donaran en gener de 2019, setembre de 2019 i gener de 2020. Necessitam calcular les pendents de les rectes d’entre setembre de 2018 i gener de 2020, per poder fer la predicció i calcular l’error.
Per calcular les pendents de les noves rectes farem la mitjana de les pendents anteriors. De les pendents ja calculades en el model obviarem la primera i la darrera, ja que no són vàlides. Per tant,
La pendent de 9-2018 a 1-2019 i de 9-2019 a 1-2020 serà : (-67.495 - 78.019)/2 = -156.115
La pendent de 1-2019 a 9-2019 serà : (61.003 + 69.522 + 57.015)/3 = 97.39
Ara, seguint el mateix procediment que en el cas de Mallorca, els nous punts d’intersecció són: 6755.746, -3384.454 i 8783.786.
Llavors l’equació per a la predicció és:
\(\hat{y}= \left\{ \begin{array}{lcc} -156.115x + 6755.746 & si & x \leq \psi_7 \\ \\ 97.39x - 3384.454 & si & \psi_7 < x \leq \psi_8 \\ \\ -156.115x + 8783.786 & si & \psi_8 < x \\ \end{array} \right.\)
on \(\psi_7 = 40\) i \(\psi_8 = 47\).
Visualitzem-ho:
p1_serie_men <- ggplot(men, aes(x=x, y=y)) + geom_line()+
ggtitle('Predicció de Menorca amb el model de regressió segmentada \nabans de la COVID') +
xlab('Índex del mes') +
ylab('Despeses mensuals en €')+
ylim(c(0,1300))
my.model1_men <- data.frame(x=1:39, y=fitted(serie1_men.seg)) #model nou amb els valors ajustats de la regressió segmentada
p1_serie_men <- p1_serie_men + geom_line(data=my.model1_men, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada
my.lines1_men <- serie1_men.seg$psi[,2]# guardam on estan els punts de canvi estimats
p1_serie_men <- p1_serie_men + geom_vline(xintercept=my.lines1_men, linetype="dashed")
p1_serie_men <- p1_serie_men + geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic
p1_serie_men <- p1_serie_men + geom_abline(intercept = 6755.746, slope=-156.115, colour="green") +
geom_abline(intercept = -3384.454, slope=97.39, colour="blue") +
geom_abline(intercept = 8783.786, slope=-156.115, colour="orange")
ggplotly(p1_serie_men)Calculem l’error de la predicció:
o1_men<-c(serie_men[40:51]) #observacions reals de la sèrie de gener de 2019 a desembre de 2019
v1_men=c(40:48)
f1_men <- sapply(v1_men, function(x) 97.39*x-3384.454) #predicció de gener 2019 a setembre 2019
v2_men=c(49:51)
f2_men <- sapply(v2_men, function(x) -156.115*x + 8783.786 )#predicció de setembre 2019 a desembre 2019
p1_men<- c(f1_men,f2_men) #predicció de 2019
sqrt(sum((p1_men-o1_men)^2)/12) #error de la predicció## [1] 203.5378- DECOMPOSE
Recordem la sèrie
plot.ts(serie1_men, main= 'Menorca abans de la COVID', ylab='Despeses mensuals en €',xlab='Any')
Ja hem dit anteriorment que té una tendència mínima i podem observar també que no hi ha variabilitat.
El mètode clàssic de descomposició separa la sèrie en subseries corresponents a la tendència, la estacionalitat i el renou.
Aquestes components es poden combinar de manera additiva o multiplicativa. En el nostre cas utilitzam el model additiu: \(y_t = \mu_t + S_t + a_t\)
decompose_s1_men<-decompose(serie1_men)
plot(decompose_s1_men, xlab="Any")
El decompose, per calcular aquestes noves tendències, el
que fa és agafar les sis tendències anteriors i les sis posteriors de la
sèrie original i en fa la mitjana. És per això que la primera que
obtenim és en abril de 2016 i la darrera en juny de 2018. Notem que per
a la predicció ens quedarem amb el valor de la darrera tendència del
decompose.
t1_men <- decompose_s1_men$trend #tendències de la sèrie sense el tros a predir abans de la COVID
t1_men## Jan Feb Mar Apr May Jun Jul Aug
## 2015
## 2016 NA NA NA 720.7354 723.6046 715.7458 707.0200 710.3858
## 2017 715.7379 724.1704 726.8771 726.2737 734.8071 744.5700 749.6183 748.9800
## 2018 753.2433 755.7308 759.9658 769.7688 777.2221 776.3413 NA NA
## Sep Oct Nov Dec
## 2015 NA NA NA
## 2016 706.4712 696.1796 698.4267 706.9596
## 2017 748.7708 749.4954 748.3637 750.5696
## 2018 NA NA NA NAEls valors de les components estacionals els calcula fent la mitjana per mesos, és a dir, per calcular la componen de gener, agafa els valors de tots els geners que tenim a la sèrie original i en fa la mitjana. Per tant, només tenim 12 valors, un per a cada mes.
s1_men <- decompose_s1_men$seasonal
s1_men <- s1_men[4:15] #estacionalitat de gener a desembre
s1_men## [1] -244.4097 -256.5797 -203.0405 -102.6617 115.6563 163.6352 243.0868
## [8] 248.9930 223.4849 188.0684 -121.3293 -254.9037Anem a fer la predicció d’aquest model
a1_men <- c(s1_men[7:12],s1_men) #estacionalitat de juliol 2018 a desembre 2019 (es per poder fer la predicció)
pred1_decompose_men <- sapply(a1_men, function(x) 776.3413 + x) #predicció de juliol de 2018 a desembre 2019p1_dec_men<-c(serie1_men[1:33], pred1_decompose_men)
A1_men<- data.frame("x" = men[1:51,]$x, "y" = men[1:51,]$y, "p"= p1_dec_men)
#ho dibuixam
grafica1_men_dec <- ggplot(A1_men)+
ylim(c(300,1100))+
geom_line(aes(x,p), color="red")+
geom_line(aes(x,y))+
geom_vline(xintercept=0) + geom_hline(yintercept=300)+
labs(title="Predicció de Menorca abans de la COVID amb el model de descomposició", x="Índex del mes", y="Despeses mensuals en €")+
theme(panel.background = element_blank())
grafica1_men_dec
Podem observar que la predicció pareix bastant bona. Calculem l’error que es comet.
sqrt(sum((c(serie_men[40:51]- pred1_decompose_men[7:18]))^2)/12) #error predicció de gener a desembre de 2019## [1] 65.76087Notem que també tenim una altra instrucció a R per fer prediccions
d’una sèrie, la funció predict(). Aquesta és basa en un
model ETS, anem a veure la predicció:
prediccio1_men <- predict(serie1_men,h=12)
plot(prediccio1_men, xlab="Any", ylab="Despeses mensuals en €")
Pareix que la predicció és bastant bona, ja que el cicle predit segueix un mateix patró que els anteriors. Anem a veure la predicció juntament amb la sèrie original:
df_prediccio1_men <- data.frame(prediccio1_men)
p1_ets_men <- data.frame("x"= 1:51, "PointForecast"=c(serie1_men,df_prediccio1_men$Point.Forecast), "Lo80" = c(rep(NA,39),df_prediccio1_men$Lo.80), "Hi80" = c(rep(NA,39),df_prediccio1_men$Hi.80), "Lo95" = c(rep(NA,39),df_prediccio1_men$Lo.95),"Hi95" = c(rep(NA,39),df_prediccio1_men$Hi.95))
grafica_pred1_ets <- ggplot((men[1:51,]))+
geom_ribbon(data = p1_ets_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data =p1_ets_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = p1_ets_men, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció del model ETS(A,N,A) a Menorca \nabans de la COVID", x="Índex del mes", y="Despeses mensuals en €") +
geom_vline(xintercept=0) + geom_hline(yintercept=0)+ theme(panel.background = element_blank())
grafica_pred1_ets
Podem observar que la predicció és bastant bona, ja que continua seguint un mateix patró. I l’error comés és d’uns 73 euros.
accuracy(prediccio1_men,serie2_men)## ME RMSE MAE MPE MAPE MASE
## Training set 7.725199 66.13662 54.76033 -0.033064 8.005419 0.7094055
## Test set -22.966948 72.64049 61.40002 -5.857499 10.400781 0.7954210
## ACF1 Theil's U
## Training set 0.05953999 NA
## Test set 0.05881637 0.5508193- SARIMA
Anem a veure quin model obtenim considerant un model estacional.
Pel que hem vist anteriorment, podem considerar que no hi ha tendència, llavors no fa falta fer cap diferència a la part regular, no obstant, sí que cal fer una diferenciació d’orde 12. Recordem que l’ACF i el PACF son:
par(mfrow=c(1,2))
acf(serie1_men)
pacf(serie1_men)
Per a la part regular obtenim un ARIMA(1,0,2). Feim una diferenciació a la part estacional, és a dir, d’ordre 12:
serie1_men_dif <- diff(serie1_men,12)
plot(serie1_men_dif, main="Sèrie sense estacionalitat", xlab="Any", ylab="Sèrie diferenciada")
Aquesta és la sèrie sense estacionalitat ni tendència, vegem com es modifiquen l’ACF i el PACF.
par(mfrow=c(2,1))
acf(serie1_men_dif, lag=36)
pacf(serie1_men_dif,lag=36)
Per a la part estacional obtenim que P=0, D=1 i Q=0.
És a dir, obtenim un ARIMA(1,0,2)(0,1,0)[12], vegem quin model ens proposa R:
auto.arima(serie1_men)## Series: serie1_men
## ARIMA(0,0,0)(0,1,0)[12]
##
## sigma^2 = 8644: log likelihood = -160.68
## AIC=323.37 AICc=323.53 BIC=324.66R ens proposa un ARIMA(0,0,0)(0,1,0)[12]. Per tant les propostes de model ARIMA són:
model1_men<-arima(serie1_men, order=c(1,0,2), seasonal = list(order=c(0,1,0), period=12))
model2_men <- arima(serie1_men, order=c(0,0,0), seasonal = list( order=c(0,1,0), period=12))Amb uns errors de 76 i 77 euros cada un.
accuracy(model1_men)## ME RMSE MAE MPE MAPE MASE
## Training set 16.36758 76.16479 51.76796 0.7966955 7.986194 0.4518997
## ACF1
## Training set -0.05174254accuracy(model2_men)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 17.76993 77.3574 53.66121 0.9457849 8.259738 0.4684264 0.07328248Anem a fer la predicció de 2019 d’aquests models i visualitzem-les.
fc_m1_men <-forecast(model1_men, h=12)
fc_m2_men <- forecast(model2_men,h=12)# model 1
pre1_men <- data.frame("Point Forecast" = serie1_men, "Lo 80" = rep(NA,39), "Hi 80"= rep(NA,39), "Lo 95" = rep(NA,39), "Hi 95" = rep(NA,39)) #dades abans de la predicció amb NA als intervals ja que només els volem per la predicció
pred_m1_men <-data.frame(fc_m1_men)
grafica_m1_men <- data.frame("x" = 1:51, "PointForecast" = c(pre1_men$Point.Forecast,pred_m1_men$Point.Forecast), "Lo80" = c(pre1_men$Lo.80, pred_m1_men$Lo.80), "Hi80"= c(pre1_men$Hi.80, pred_m1_men$Hi.80), "Lo95" = c(pre1_men$Lo.95, pred_m1_men$Lo.95), "Hi95" = c(pre1_men$Hi.95, pred_m1_men$Hi.95))
grafica1_men <- ggplot(men[1:51,])+
geom_ribbon(data = grafica_m1_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m1_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m1_men, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció de Menorca amb el model ARIMA(1,0,2)(0,1,0)[12] \nabans de la COVID", x="Índex del mes", y="Despeses mensuals en €") +
geom_vline(xintercept=0) + geom_hline(yintercept=0)+ theme(panel.background = element_blank())
grafica1_men
#model 2
pred_m2_men <-data.frame(fc_m2_men)
grafica_m2_men <- data.frame("x" = 1:51, "PointForecast" = c(pre1_men$Point.Forecast,pred_m2_men$Point.Forecast), "Lo80" = c(pre1_men$Lo.80, pred_m2_men$Lo.80), "Hi80"= c(pre1_men$Hi.80, pred_m2_men$Hi.80), "Lo95" = c(pre1_men$Lo.95, pred_m2_men$Lo.95), "Hi95" = c(pre1_men$Hi.95, pred_m2_men$Hi.95))
grafica2_men <- ggplot(men[1:51,])+
geom_ribbon(data = grafica_m2_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m2_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m2_men, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció de Menorca amb el model ARIMA (0,0,0)(0,1,0)[12] \nabans de la COVID", x="Índex del mes", y="Despeses mensuals en €") +
geom_vline(xintercept=0) + geom_hline(yintercept=0)+ theme(panel.background = element_blank())
grafica2_men
Vegem i estudiem els errors de la predicció:
accuracy(fc_m1_men, serie_men[40:51], h=12)## ME RMSE MAE MPE MAPE MASE
## Training set 16.36758 76.16479 51.76796 0.7966955 7.986194 0.4518997
## Test set -18.96009 80.00195 54.61251 -3.5716422 9.554005 0.4767307
## ACF1
## Training set -0.05174254
## Test set NAaccuracy(fc_m2_men,serie_men[40:51], h=12)## ME RMSE MAE MPE MAPE MASE
## Training set 17.76993 77.3574 53.66121 0.9457849 8.259738 0.4684264
## Test set -19.24750 80.4762 55.42417 -3.6350339 9.727838 0.4838159
## ACF1
## Training set 0.07328248
## Test set NAcheckresiduals(fc_m1_men)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,2)(0,1,0)[12]
## Q* = 5.4139, df = 5, p-value = 0.3675
##
## Model df: 3. Total lags used: 8checkresiduals(fc_m2_men)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,0,0)(0,1,0)[12]
## Q* = 6.1129, df = 8, p-value = 0.6346
##
## Model df: 0. Total lags used: 8par(mfrow=c(1,2))
qqPlot(fc_m1_men$residuals, id=FALSE, mean=mean(fc_m1_men$residuals), sd=sd(fc_m1_men$residuals))
qqPlot(fc_m2_men$residuals, id=FALSE, mean=mean(fc_m2_men$residuals), sd=sd(fc_m2_men$residuals))
shapiro.test(fc_m1_men$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m1_men$residuals
## W = 0.91595, p-value = 0.006534shapiro.test(fc_m2_men$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m2_men$residuals
## W = 0.93545, p-value = 0.02697lillie.test(fc_m1_men$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m1_men$residuals
## D = 0.21097, p-value = 0.0001394lillie.test(fc_m2_men$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m2_men$residuals
## D = 0.17916, p-value = 0.002857Resum dels errors que cometen cada un dels models en aquest cas:
| reg. segmentada | descomposició clàssica | ETS(A,N,A) | ARIMA (1,0,2)(0,1,0)[12] | ARIMA (0,0,0)(0,1,0)[12] | |
|---|---|---|---|---|---|
| error model | 70.71 | 66.14 | 76.16 | 77.36 | |
| error predicció | 203.76 | 65.76 | 72.64 | 80 | 80.48 |
2.2. DURANT LA COVID
- REGRESSIÓ SEGMENTADA
El valor d’\(R^2\) per a la regressió lineal és:
serie2_men.lm <- lm(y~x, data= men[1:51,])
summary(serie2_men.lm)$adj.r.squared## [1] 0.01782343Aleshores utilitzem el paquet segmented per ajustar les
nostres dades a una regressió lineal segmentada i millorar aquest
valor.
Vegem els punts de canvi:
punts_canvi_serie2_men <-selgmented(serie2_men.lm, Kmax = 10, type="bic")## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 .. 7 .. 8 .. 9 .. 10 ..
##
## BIC to detect no. of breakpoints:
## 0 1 2 3 4 5 6 7
## 557.4237 564.4217 568.4880 572.6311 580.3316 571.0819 573.0862 564.6621
## 8 9 10
## 565.2093 524.6382 530.1464
##
## No. of selected breakpoints: 8 (1 breakpoint(s) removed due to small slope diff)serie2_men.seg <- segmented(serie2_men.lm, seg.Z=~x, psi = c(5,10,14,21,29,35,40,47))
summary(serie2_men.seg)$psi## Initial Est. St.Err
## psi1.x 5 3.954453 0.9076704
## psi2.x 10 11.839275 0.4112265
## psi3.x 14 15.390922 0.3938302
## psi4.x 21 22.697828 0.4783673
## psi5.x 29 28.566720 0.4625795
## psi6.x 35 35.431934 0.3948826
## psi7.x 40 40.285772 0.3812246
## psi8.x 47 47.349153 0.4179509Aquests són en: 1-2016, 9-2016, 12-2016, 8-2017, 2-2018, 8-2018, 1-2019 i 8-2019.
Ara el valor de \(R^2\) de la segmentada és
summary(serie2_men.seg)$adj.r.squared## [1] 0.8254133Anem a visualitzar la regressió segmentada damunt les nostres dades
p_serie2_men <- ggplot(men[1:51,], aes(x=x, y=y)) + geom_line()+
ylim(c(0,1300)) +
ggtitle('Regressió lineal i segmentada sobre les dades \nde Menorca durant la COVID') +
xlab('índex del mes')+
ylab('Despeses mensuals en €')
my.coef2_men <- coef(serie2_men.lm) #coeficients de la recta de regressió lineal
p_serie2_men <- p_serie2_men + geom_abline(intercept=my.coef2_men[1], slope=my.coef2_men[2], color="green") #afegim la recta de regressió lineal
my.model2_men <- data.frame(x=1:51, y=fitted(serie2_men.seg)) #model nou amb els valors ajustats de la regressió segmentada
p_serie2_men <- p_serie2_men + geom_line(data=my.model2_men, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada
my.lines2_men <- serie2_men.seg$psi[,2]# guardam on estan els punts de canvi estimats
p_serie2_men <- p_serie2_men + geom_vline(xintercept=my.lines2_men, linetype="dashed")
p_serie2_men <- p_serie2_men + geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
ggplotly(p_serie2_men)Per poder escriure la funció necessitam les pendents i els punts d’intersecció amb l’eix OY, que ens ho donen les següents funcions:
#PUNTS D'INTERSECCIÓ
intercept(serie2_men.seg)## $x
## Est.
## intercept1 841.24
## intercept2 171.78
## intercept3 3347.90
## intercept4 -1263.70
## intercept5 3720.40
## intercept6 -2821.90
## intercept7 7248.30
## intercept8 -4337.70
## intercept9 9635.50#PENDENTS
slope(serie2_men.seg)## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -97.165 66.526 -1.4605 -232.510 38.184
## slope2 72.128 14.517 4.9684 42.593 101.660
## slope3 -196.140 42.075 -4.6616 -281.740 -110.540
## slope4 103.490 17.780 5.8208 67.320 139.670
## slope5 -116.090 22.490 -5.1620 -161.850 -70.337
## slope6 112.920 17.780 6.3512 76.750 149.100
## slope7 -171.290 29.751 -5.7573 -231.820 -110.760
## slope8 116.310 17.780 6.5415 80.133 152.480
## slope9 -178.800 42.075 -4.2496 -264.400 -93.200L’error del model de regressió segmentada és (ens interessa el RMSE):
accuracy(serie2_men.seg)## ME RMSE MAE MPE MAPE MASE
## Training set 0 75.67986 63.29412 -1.310022 9.753939 0.3068017Anem a fer la predicció del 2020:
Recordem punts de canvi en: 1-2016, 9-2016, 12-2016, 8-2017, 2-2018, 8-2018, 1-2019, 8-2019. Aleshores, seguint el mateix criteri que abans, tenim que les noves previsions seran en 1-2020, 8-2020 i 1-2021. De la mateixa manera que a Mallorca treim les noves pendents i punts d’intercessció i obtenim que la funció pel tros predit és:
\(\hat{y}= \left\{ \begin{array}{lcc} -161.173x + 8801.12 & si & x \leq \psi_9 \\ \\ 101.212x - 4842.9 & si & \psi_9 < x \leq \psi_{10} \\ \\ -161.173x + 10637.82 & si & \psi_{10} < x \\ \end{array} \right.\)
on \(\psi_9 = 52\) i \(\psi_8 = 59\).
Visualitzem-ho:
#ho graficam
p2_serie_men <- ggplot(men, aes(x=x, y=y)) + geom_line()+
ggtitle('Predicció de Menorca amb el model \nde regressió segmentada durant la COVID') +
xlab('Índex del mes') +
ylab('Despeses mensuals en €') +
ylim(c(0,1300))
my.model2_men <- data.frame(x=1:51, y=fitted(serie2_men.seg)) #model nou amb els valors ajustats de la regressió segmentada
p2_serie_men <- p2_serie_men + geom_line(data=my.model2_men, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada
my.lines2_men <- serie2_men.seg$psi[,2]# guardam on estan els punts de canvi estimats
p2_serie_men <- p2_serie_men + geom_vline(xintercept=my.lines2_men, linetype="dashed") #línies verticals en els punts de canvi
p2_serie_men <- p2_serie_men + geom_vline(xintercept=0) + geom_hline(yintercept=0)+ theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic
p2_serie_men <- p2_serie_men + geom_abline(intercept = 8801.12, slope=-161.173, colour="green") +
geom_abline(intercept = -4842.9 , slope=101.212, colour="blue") +
geom_abline(intercept = 10637.82 , slope=-161.173, colour="orange")
ggplotly(p2_serie_men)Calculem l’error de la predicció:
o2_men<-c(serie_men[52:63]) # dades reals per fer predicció del 2020
v3_men=c(52:59)
f3_men <- sapply(v3_men, function(x) 101.212*x-4842.9) #predicció de gener 2020 a agost 2020
v4_men=c(60:63)
f4_men <- sapply(v4_men, function(x) -161.173*x + 10637.82) #predicció de setembre 2020 a desembre 2020
p2_men <- c(f3_men,f4_men) #predicció de 2020
sqrt(sum((p2_men-o2_men)^2)/12)## [1] 354.9417- DECOMPOSE
La descomposició de la sèrie en aquest cas és la següent
decompose_s2_men <- decompose(serie2_men)
plot(decompose_s2_men, xlab="Any")
On les components de tendència són:
t2_men <- decompose_s2_men$trend #tendències de la sèrie sense el tros a predir abans de la COVID
t2_men## Jan Feb Mar Apr May Jun Jul Aug
## 2015
## 2016 NA NA NA 720.7354 723.6046 715.7458 707.0200 710.3858
## 2017 715.7379 724.1704 726.8771 726.2737 734.8071 744.5700 749.6183 748.9800
## 2018 753.2433 755.7308 759.9658 769.7688 777.2221 776.3413 772.0879 771.4596
## 2019 762.4079 762.7075 763.6933 761.6817 750.9704 749.0288 NA NA
## Sep Oct Nov Dec
## 2015 NA NA NA
## 2016 706.4712 696.1796 698.4267 706.9596
## 2017 748.7708 749.4954 748.3637 750.5696
## 2018 768.9025 765.1746 763.9954 762.3817
## 2019 NA NA NA NAI les estacionals:
s2_men <- decompose_s2_men$seasonal #estacionalitat
s2_men <- s2_men[4:15] #estacionalitat de gener a desembre
s2_men## [1] -262.68892 -267.21878 -232.48462 -120.92743 115.23893 173.50351
## [7] 256.53538 256.21233 228.83927 226.40427 -97.99781 -275.41614Fem la predicció:
a2_men <- c(s2_men[7:12],s2_men) #estacionalitat de juliol 2019 a desembre 2029 (es per poder fer la predicció)
pred2_decompose_men <- sapply(a2_men, function(x) 749.0288 + x) #predicció de juliol de 2019 a desembre 2020
p2_dec_men<-c(serie2_men[1:45], pred2_decompose_men)
A2_men<- data.frame("x" = men[1:63,]$x, "y" = men[1:63,]$y, "p"= p2_dec_men)
grafica2_men_dec <- ggplot(A2_men)+
geom_line(aes(x,p), color="red")+
geom_line(aes(x,y))+
labs(title="Predicció de Menorca durant la COVID amb el model de descomposició", x="Índex del mes", y="Despesese mensuals en €")+
geom_vline(xintercept = 0)+ geom_hline(yintercept = 0)+ theme(panel.background = element_blank())
grafica2_men_dec
L’error que es comet és:
sqrt(sum((c(serie_men[52:63]- pred2_decompose_men[7:18]))^2)/12)## [1] 346.3806Vegem, igual que abans, la predicció amb la funció
predict():
prediccio2_men <- predict(serie2_men,h=12)
plot(prediccio2_men, xlab="Any", ylab ="Despeses menusals en €")
Vegem com queda la predicció sobre la sèrie original:
df_prediccio2_men <- data.frame(prediccio2_men)
p2_ets_men <- data.frame("x"= 1:63, "PointForecast"=c(serie2_men,df_prediccio2_men$Point.Forecast), "Lo80" = c(rep(NA,51),df_prediccio2_men$Lo.80), "Hi80" = c(rep(NA,51),df_prediccio2_men$Hi.80), "Lo95" = c(rep(NA,51),df_prediccio2_men$Lo.95),"Hi95" = c(rep(NA,51),df_prediccio2_men$Hi.95))
grafica_pred2_ets <- ggplot((men[1:63,]))+
geom_ribbon(data = p2_ets_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data =p2_ets_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = p2_ets_men, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció del model ETS(A,N,A) a Menorca durant la COVID", y="Despeses mensuals en €", x="Índex del mes") +
geom_vline(xintercept = 0) + geom_hline(yintercept = 0)+ theme(panel.background = element_blank())
grafica_pred2_ets
Si calculam l’error comés, aquest és d’uns 351 euros.
accuracy(prediccio2_men, serie3_men)## ME RMSE MAE MPE MAPE MASE
## Training set 10.85421 64.23595 53.31916 0.5545536 7.951054 0.7563633
## Test set -213.59590 350.65163 240.79195 -Inf Inf 3.4157744
## ACF1 Theil's U
## Training set 0.1249660 NA
## Test set 0.4522522 0- SARIMA
Quin model proposam nosaltres? Vegem l’ACF i el PACF.
par(mfrow=c(1,2))
acf(serie2_men)
pacf(serie2_men)
Per a la part regular obtenim un ARIMA(4,0,2).
Observam que hi ha estacionalitat, llavors feim una diferenciació d’ordre 12.
serie2_men_diff <- diff(serie2_men,12)
plot(serie2_men_diff, main="Sèrie sense estacionalitat", xlab="Any", ylab="Sèrie diferenciada")
Ara ja no s’observa l’estacionalitat, llavors hem de fer una diferenciació D=1. Vegem quins són els nous ACF i PACF.
par(mfrow=c(1,2))
acf(serie2_men_diff,lag=36)
pacf(serie2_men_diff,lag=36)
Obtenim que P=0 i Q=0. Per tant el model que nosaltres proposam és un ARIMA(4,0,2)(0,1,0)[12].
R proposa el següent model:
auto.arima(serie2_men)## Series: serie2_men
## ARIMA(0,0,0)(0,1,0)[12]
##
## sigma^2 = 7977: log likelihood = -230.53
## AIC=463.06 AICc=463.17 BIC=464.73Llavors els models que tenim són:
model3_men <- arima(serie2_men, order=c(4,0,2), seasonal = list(order=c(0,1,0), period=12))
model4_men <- arima(serie2_men, order=c(0,0,0), seasonal = list( order=c(0,1,0), period=12))I els seus errors:
accuracy(model3_men)## ME RMSE MAE MPE MAPE MASE
## Training set 8.802988 69.36035 46.9251 0.06608693 7.264296 0.4210999
## ACF1
## Training set -0.007604448accuracy(model4_men)## ME RMSE MAE MPE MAPE MASE
## Training set 9.059943 78.10244 54.07602 -0.1320548 8.605173 0.4852714
## ACF1
## Training set 0.02022333Vegem les prediccions
fc_m3_men <- forecast(model3_men, h=12)
fc_m4_men <- forecast(model4_men,h=12)#primer model
pre2_men <- data.frame("Point Forecast" = serie2_men, "Lo 80" = rep(NA,51), "Hi 80"= rep(NA,51), "Lo 95" = rep(NA,51), "Hi 95" = rep(NA,51))
pred_m3_men <-data.frame(fc_m3_men)
grafica_m3_men <- data.frame("x" = 1:63, "PointForecast" = c(pre2_men$Point.Forecast,pred_m3_men$Point.Forecast), "Lo80" = c(pre2_men$Lo.80, pred_m3_men$Lo.80), "Hi80"= c(pre2_men$Hi.80, pred_m3_men$Hi.80), "Lo95" = c(pre2_men$Lo.95, pred_m3_men$Lo.95), "Hi95" = c(pre2_men$Hi.95, pred_m3_men$Hi.95))
grafica3_men <- ggplot(men[1:63,])+
geom_ribbon(data = grafica_m3_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m3_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m3_men, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Prediccó de Menorca amb el model ARIMA (4,0,2)(0,1,0)[12] \ndurant la COVID", x="Índex del mes", y="Despeses mensuals en €")+
geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica3_men
#segon model
pred_m4_men <-data.frame(fc_m4_men)
grafica_m4_men <- data.frame("x" = 1:63, "PointForecast" = c(pre2_men$Point.Forecast,pred_m4_men$Point.Forecast), "Lo80" = c(pre2_men$Lo.80, pred_m4_men$Lo.80), "Hi80"= c(pre2_men$Hi.80, pred_m4_men$Hi.80), "Lo95" = c(pre2_men$Lo.95, pred_m4_men$Lo.95), "Hi95" = c(pre2_men$Hi.95, pred_m4_men$Hi.95))
grafica4_men <- ggplot(men[1:63,])+
geom_ribbon(data = grafica_m4_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m4_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m4_men, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Prediccó de Menorca amb el model ARIMA (0,0,0)(0,1,0)[12] \ndurant la COVID", x="Índex del mes", y="Despeses mensuals en €")+
geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica4_men
Vegem i estudiem els seus errors:
accuracy(fc_m3_men,serie_men[52:63], h=12)## ME RMSE MAE MPE MAPE MASE
## Training set 8.802988 69.36035 46.9251 0.06608693 7.264296 0.4210999
## Test set -207.871895 347.65796 256.7831 -Inf Inf 2.3043392
## ACF1
## Training set -0.007604448
## Test set NAaccuracy(fc_m4_men,serie_men[52:63], h=12)## ME RMSE MAE MPE MAPE MASE
## Training set 9.059943 78.10244 54.07602 -0.1320548 8.605173 0.4852714
## Test set -205.551667 350.77019 256.22500 -Inf Inf 2.2993307
## ACF1
## Training set 0.02022333
## Test set NAcheckresiduals(fc_m3_men)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(4,0,2)(0,1,0)[12]
## Q* = 1.9353, df = 4, p-value = 0.7477
##
## Model df: 6. Total lags used: 10checkresiduals(fc_m4_men)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,0,0)(0,1,0)[12]
## Q* = 4.2132, df = 10, p-value = 0.9372
##
## Model df: 0. Total lags used: 10par(mfrow=c(1,2))
qqPlot(fc_m3_men$residuals, id=FALSE, mean=mean(fc_m3_men$residuals), sd=sd(fc_m3_men$residuals))
qqPlot(fc_m4_men$residuals, id=FALSE, mean=mean(fc_m3_men$residuals), sd=sd(fc_m3_men$residuals))
shapiro.test(fc_m3_men$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m3_men$residuals
## W = 0.94617, p-value = 0.02186shapiro.test(fc_m4_men$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m4_men$residuals
## W = 0.9436, p-value = 0.01717lillie.test(fc_m3_men$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m3_men$residuals
## D = 0.16187, p-value = 0.001878lillie.test(fc_m4_men$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m4_men$residuals
## D = 0.17622, p-value = 0.0004126Resum dels errors que comet cada un del models anteriors:
| reg. segmentada | descomposició clàssica | ETS(A,N,A) | ARIMA (4,0,2)(0,1,0)[12] | ARIMA (0,0,0)(0,1,0)[12] | |
|---|---|---|---|---|---|
| error model | 75.68 | 64.24 | 69.36 | 78.1 | |
| error predicció | 354.9417 | 346.38 | 350.65 | 347.66 | 350.77 |
2.3. DESPRÉS DE LA COVID
- REGRESSIÓ SEGMENTADA
Anem a fer ús del paquet segmented.
Vegem els punts de canvi que, igual que a Mallorca, en aquest cas ens indica que no n’hi ha:
serie3_men.lm <- lm(y~x, data=men[1:63,])
punts_canvi_serie3_men <-selgmented(serie3_men.lm, Kmax = 10, type="bic")## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 ..
## (search truncated at 6 breakpoints due to increasing values of BIC)
##
## BIC to detect no. of breakpoints:
## 0 1 2 <NA> 3 4 5
## 700.6072 704.9516 712.8201 713.8201 714.3654 720.6737 720.9406
##
## No. of selected breakpoints: 0serie3_men.seg <- segmented(serie3_men.lm, seg.Z=~x, psi = c(5,10,16,22,29,36,39,45,55,58))
summary(serie3_men.seg)$psi## Initial Est. St.Err
## psi1.x 5 4.206349 0.7961087
## psi2.x 10 11.543480 0.5612034
## psi3.x 16 16.046980 0.4894338
## psi4.x 22 22.578797 0.5367818
## psi5.x 29 28.338906 0.5479134
## psi6.x 36 36.541203 0.3634258
## psi7.x 39 39.431366 0.3513731
## psi8.x 45 46.731891 0.4875429
## psi9.x 55 55.984004 0.2442491
## psi10.x 58 58.007390 0.3267941Ens queden els punts de canvi a: 1-2016, 9-2016, 1-2017, 8-2017, 1-2018, 10-2018, 12-2018, 8-2019, 5-2020 i 7-2020.
Ara, el valor de \(R^2\) de la segmentada és
summary(serie3_men.seg)$adj.r.squared## [1] 0.813196Anem a visualitzar la regressió segmentada damunt les nostres dades
p_serie3_men <- ggplot(men[1:63,], aes(x=x, y=y)) + geom_line()+
ylim(c(0,1300))+
ggtitle('Regressió lineal i segmentada sobre les dades de \nMenorca després de la COVID') +
xlab('índex del mes')+
ylab('Despeses mensuals en €')
my.coef3_men <- coef(serie3_men.lm) #coeficients de la recta de regressió lineal
p_serie3_men <- p_serie3_men + geom_abline(intercept=my.coef3_men[1], slope=my.coef3_men[2], color="green") #afegim la recta de regressió lineal
my.model3_men <- data.frame(x=1:63, y=fitted(serie3_men.seg)) #model nou amb els valors ajustats de la regressió segmentada
p_serie3_men <- p_serie3_men + geom_line(data=my.model3_men, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada
my.lines3_men <- serie3_men.seg$psi[,2]# guardam on estan els punts de canvi estimats
p_serie3_men <- p_serie3_men + geom_vline(xintercept=my.lines3_men, linetype="dashed")
p_serie3_men <- p_serie3_men + geom_vline(xintercept = 0)+ geom_hline(yintercept = 0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic
ggplotly(p_serie3_men)Per escriure la funció a trossos necessitam les pendents i els punts d’intersecció de les rectes:
#punts d'intersecció
intercept(serie3_men.seg)## $x
## Est.
## intercept1 824.68
## intercept2 139.10
## intercept3 2722.00
## intercept4 -1500.40
## intercept5 3720.40
## intercept6 -2214.10
## intercept7 12566.00
## intercept8 -4119.80
## intercept9 6348.40
## intercept10 -22225.00
## intercept11 4038.00#pendents
slope(serie3_men.seg)## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -87.227 47.552 -1.8344 -183.260 8.8053
## slope2 75.759 20.094 3.7702 35.178 116.3400
## slope3 -148.000 33.624 -4.4016 -215.900 -80.0930
## slope4 115.130 25.417 4.5297 63.802 166.4600
## slope5 -116.090 25.417 -4.5675 -167.420 -64.7620
## slope6 93.318 16.407 5.6877 60.183 126.4500
## slope7 -311.160 75.186 -4.1385 -463.000 -159.3100
## slope8 112.000 20.094 5.5738 71.421 152.5800
## slope9 -112.000 13.727 -8.1594 -139.730 -84.2810
## slope10 398.380 75.186 5.2987 246.540 550.2300
## slope11 -54.371 33.624 -1.6170 -122.280 13.5340L’error de la regressió segmentada és:
accuracy(serie3_men.seg)## ME RMSE MAE MPE MAPE MASE
## Training set 1.352864e-15 85.77708 67.2011 -Inf Inf 0.3253227Anem a fer la predicció per a l’any 2021. Recordem que els punts de canvi es donen a 1-2016, 9-2016, 1-2017, 8-2017, 1-2018, 10-2018, 12-2018, 8-2019, 5-2020 i 7-2020. Degut a la pertorbació del COVID, sí que hi segueix havent un punt de canvi en estiu i l’altre en hivern successivament però ara no és donen al mateix mes. Per això, per predir els següents punts de canvi agafarem el mes més freqüent. Aleshores els següents punts de canvi seran en 1-2021, 8-2021 I 1-2022.
- La pendent de 7-2020 a 1-2021 i de 8-2021 a 1-2022 serà : -171.8125
- La pendent de 1-2021 a 8-2021 serà : 158.9174
Els nous punts d’intersecció es calculen de forma anàloga que a Mallorca. Per a les tres noves rectes obtenim que són \(n_1\)=10851.87 , \(n_2\)=−10314.85 i \(n_3\)=13166.98.
Aleshores la predicció de 2021 serà:
\(\hat{y}= \left\{ \begin{array}{lcc} -171.813x + 10851.87 & si & x \leq \psi_{11} \\ \\ 158.9174x - 10314.85 & si & \psi_{11} < x \leq \psi_{12} \\ \\ -171.8125x + 13166.98 & si & \psi_{12} < x \\ \end{array} \right.\)
on \(\psi_9 = 64\) i \(\psi_8 = 71\).
Visualitzem-la:
p3_serie_men <- ggplot(men, aes(x=x, y=y)) + geom_line()+
ylim(c(-200,1300)) + ggtitle('Predicció de Menorca amb el model de \nregressió segmentada després de la COVID') + xlab('índex del mes')+ ylab('Despeses mensuals en €')
my.model3_men <- data.frame(x=1:63, y=fitted(serie3_men.seg)) #model nou amb els valors ajustats de la regressió segmentada
p3_serie_men <- p3_serie_men + geom_line(data=my.model3_men, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada
my.lines3_men <- serie3_men.seg$psi[,2]# guardam on estan els punts de canvi estimats
p3_serie_men <- p3_serie_men + geom_vline(xintercept=my.lines3_men, linetype="dashed")
p3_serie_men <- p3_serie_men + geom_vline(xintercept = 0)+ geom_hline(yintercept = 0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic
p3_serie_men <- p3_serie_men + geom_abline(intercept = 10851.87, slope=-171.8125, colour="green") +
geom_abline(intercept = -10314.85 , slope=158.9174, colour="blue") +
geom_abline(intercept = 13166.98 , slope=-171.8125, colour="orange")
ggplotly(p3_serie_men)Vegem quin és aquest error que es comet:
o3_men<-c(serie_men[52:63]) # dades reals per fer predicció del 2021
v5_men=c(52:59)
f5_men <- sapply(v5_men, function(x) 158.9174*x-10314.85) #predicció de gener 2021 a agost 2021
v6_men=c(60:63)
f6_men <- sapply(v6_men, function(x) -171.8125*x + 13166.98) #predicció de setembre 2021 a desembre 2021
p3_men <- c(f5_men,f6_men) #predicció de 2021
sqrt(sum((p3_men-o3_men)^2)/12)## [1] 1977.984- DECOMPOSE
Visualitzem la sèrie descomposada:
decompose_s3_men <- decompose(serie3_men)
plot(decompose_s3_men, xlab="Any")
Les components de tendència són:
t3_men <- decompose_s3_men$trend #tendències de la sèrie sense el tros a predir abans de la COVID
t3_men## Jan Feb Mar Apr May Jun Jul Aug
## 2015
## 2016 NA NA NA 720.7354 723.6046 715.7458 707.0200 710.3858
## 2017 715.7379 724.1704 726.8771 726.2737 734.8071 744.5700 749.6183 748.9800
## 2018 753.2433 755.7308 759.9658 769.7688 777.2221 776.3413 772.0879 771.4596
## 2019 762.4079 762.7075 763.6933 761.6817 750.9704 749.0288 757.8138 761.0633
## 2020 608.6500 589.2296 571.2750 548.4712 539.5446 546.7817 NA NA
## Sep Oct Nov Dec
## 2015 NA NA NA
## 2016 706.4712 696.1796 698.4267 706.9596
## 2017 748.7708 749.4954 748.3637 750.5696
## 2018 768.9025 765.1746 763.9954 762.3817
## 2019 763.7283 741.4167 680.8050 632.0783
## 2020 NA NA NA NAI les d’estacionalitat:
s3_men <- decompose_s3_men$seasonal #estacionalitat
s3_men <- s3_men[4:15] #estacionalitat de gener a desembreAnem a fer la predicció
a3_men <- c(s3_men[7:12],s3_men) #estacionalitat de juliol 2020 a desembre 2021 (es per poder fer la predicció)
pred3_decompose_men <- sapply(a3_men, function(x) 546.7817 + x) #predicció de juliol de 2019 a desembre 2020
p3_dec_men<-c(serie3_men[1:57], pred3_decompose_men)
A3_men<- data.frame("x" = men[1:75,]$x, "y" = men[1:75,]$y, "p"= p3_dec_men)
grafica3_men_dec <- ggplot(A3_men)+
geom_line(aes(x,p), color="red")+
geom_line(aes(x,y))+
labs(title="Predicció de Menorca després de la COVID amb el model de descomposició", x="Índex del mes", y="Despeses mensualns en €")+
geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica3_men_dec
Calculem l’error de la predicció
sqrt(sum((c(serie_men[64:75]- pred3_decompose_men[7:18]))^2)/12)## [1] 218.4594Amb la funció predict() la predicció seria la
següent:
prediccio3_men <- predict(serie3_men,h=12)
plot(prediccio3_men, xlab="Any", ylab="Despeses mensuals en €")
Vegem-la juntament amb les nostres dades:
df_prediccio3_men <- data.frame(prediccio3_men)
p3_ets_men <- data.frame("x"= 1:75, "PointForecast"=c(serie3_men,df_prediccio3_men$Point.Forecast), "Lo80" = c(rep(NA,63),df_prediccio3_men$Lo.80), "Hi80" = c(rep(NA,63),df_prediccio3_men$Hi.80), "Lo95" = c(rep(NA,63),df_prediccio3_men$Lo.95),"Hi95" = c(rep(NA,63),df_prediccio3_men$Hi.95))
grafica_pred3_ets <- ggplot((men[1:75,]))+
geom_ribbon(data = p3_ets_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data =p3_ets_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = p3_ets_men, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció del model ETS(A,N,A) a Menorca després de la COVID", x="Índex del mes", y="Despeses mensuals en €")+
geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica_pred3_ets
Obtenim un error d’uns 117 euros.
accuracy(prediccio3_men,serie_men)## ME RMSE MAE MPE MAPE MASE
## Training set 1.544179 126.6657 89.34926 -Inf Inf 0.7824237
## Test set -65.427782 117.3246 106.49793 -8.062414 13.97902 0.9325931
## ACF1 Theil's U
## Training set 0.06833943 NA
## Test set 0.46102653 1.084369- SARIMA
Quin model proposam nosaltres? Vegem l’ACF i el PACF:
par(mfrow=c(1,2))
acf(serie3_men)
pacf(serie3_men)
Per a la part regular obtenim p=2, q=2 i d=0.
Sí hi ha indicació d’estacionalitat, llavors feim una diferenciació d’ordre 12 i tornam a calcular l’ACF i el PACF:
serie3_men_diff <- diff(serie3_men,12)
plot.ts(serie3_men_diff, main="Sèrie sense estacionalitat", ylab="Sèrie diferenciada", xlab="Any")
par(mfrow=c(1,2))
acf(serie3_men_diff)
pacf(serie3_men_diff)
Aleshores obtenim un model ARIMA(2,0,2)(1,1,1)[12].
R proposa el següent model:
auto.arima(serie3_men)## Series: serie3_men
## ARIMA(0,1,2)(0,1,1)[12]
##
## Coefficients:
## ma1 ma2 sma1
## -0.2359 -0.5527 -0.4631
## s.e. 0.1334 0.1258 0.2793
##
## sigma^2 = 22849: log likelihood = -322.38
## AIC=652.76 AICc=653.65 BIC=660.41Llavors tenim els següents models
model5_men <- arima(serie3_men, order=c(2,0,2), seasonal = list( order=c(1,1,1), period=12))
model6_men <- arima(serie3_men, order=c(0,1,2), seasonal = list( order=c(0,1,1), period=12))Amb uns errors de:
accuracy(model5_men)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -8.492597 122.6242 77.91071 -Inf Inf 0.6601132 0.003622309accuracy(model6_men)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -12.09391 130.56 82.80016 -Inf Inf 0.70154 -0.004851061Vegem les prediccions
fc_m5_men <- forecast(model5_men, h=12)
fc_m6_men <- forecast(model6_men,h=12)#grafiquem-les
#primer model
pre3_men <- data.frame("Point Forecast" = serie3_men, "Lo 80" = rep(NA,63), "Hi 80"= rep(NA,63), "Lo 95" = rep(NA,63), "Hi 95" = rep(NA,63))
pred_m5_men <-data.frame(fc_m5_men)
grafica_m5_men <- data.frame("x" = 1:75, "PointForecast" = c(pre3_men$Point.Forecast,pred_m5_men$Point.Forecast), "Lo80" = c(pre3_men$Lo.80, pred_m5_men$Lo.80), "Hi80"= c(pre3_men$Hi.80, pred_m5_men$Hi.80), "Lo95" = c(pre3_men$Lo.95, pred_m5_men$Lo.95), "Hi95" = c(pre3_men$Hi.95, pred_m5_men$Hi.95))
grafica5_men <- ggplot(men[1:75,])+
geom_ribbon(data = grafica_m5_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m5_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m5_men, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció de Menorca amb el model ARIMA (2,0,2)(1,1,1)[12] \ndesprés de la COVID", x="Índex del mes", y="Despeses mensuals en €")+
geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica5_men
#segon model
pred_m6_men <-data.frame(fc_m6_men)
grafica_m6_men <- data.frame("x" = 1:75, "PointForecast" = c(pre3_men$Point.Forecast,pred_m6_men$Point.Forecast), "Lo80" = c(pre3_men$Lo.80, pred_m6_men$Lo.80), "Hi80"= c(pre3_men$Hi.80, pred_m6_men$Hi.80), "Lo95" = c(pre3_men$Lo.95, pred_m6_men$Lo.95), "Hi95" = c(pre3_men$Hi.95, pred_m6_men$Hi.95))
grafica6_men <- ggplot(men[1:75,])+
geom_ribbon(data = grafica_m6_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
geom_ribbon(data = grafica_m6_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
geom_line(data = grafica_m6_men, aes(x, PointForecast), colour = "blue") +
geom_line(aes(x,y), color="red")+
labs(title="Predicció de Menorca amb el model ARIMA (0,1,2)(0,1,1)[12] \ndesprés de la COVID", x="Índex del mes", y="Despeses mensuals en €")+
geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())
grafica6_men
Vegem quin és l’error de cada model i estudiem-los
accuracy(fc_m5_men,serie_men[64:75], h=12)## ME RMSE MAE MPE MAPE MASE
## Training set -8.492597 122.6242 77.91071 -Inf Inf 0.6601132
## Test set 76.160973 137.6361 89.64163 11.16486 12.81547 0.7595056
## ACF1
## Training set 0.003622309
## Test set NAaccuracy(fc_m6_men,serie_men[64:75], h=12)## ME RMSE MAE MPE MAPE MASE
## Training set -12.09391 130.5600 82.80016 -Inf Inf 0.701540
## Test set 218.06171 257.1519 218.06171 30.86593 30.86593 1.847569
## ACF1
## Training set -0.004851061
## Test set NAcheckresiduals(fc_m5_men)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,0,2)(1,1,1)[12]
## Q* = 9.1279, df = 7, p-value = 0.2436
##
## Model df: 6. Total lags used: 13checkresiduals(fc_m6_men)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,2)(0,1,1)[12]
## Q* = 12.831, df = 10, p-value = 0.2333
##
## Model df: 3. Total lags used: 13par(mfrow=c(1,2))
qqPlot(fc_m5_men$residuals, id=FALSE, mean=mean(fc_m5_men$residuals), sd=sd(fc_m5_men$residuals))
qqPlot(fc_m6_men$residuals, id=FALSE, mean=mean(fc_m6_men$residuals), sd=sd(fc_m6_men$residuals))
shapiro.test(fc_m5_men$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m5_men$residuals
## W = 0.83562, p-value = 7.234e-07shapiro.test(fc_m6_men$residuals)##
## Shapiro-Wilk normality test
##
## data: fc_m6_men$residuals
## W = 0.86346, p-value = 4.995e-06lillie.test(fc_m5_men$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m5_men$residuals
## D = 0.21125, p-value = 2.001e-07lillie.test(fc_m6_men$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: fc_m6_men$residuals
## D = 0.15808, p-value = 0.0004721Resum dels errors que comet cada un del models anteriors:
| reg. segmentada | descomposició clàssica | ETS(A,N,A) | ARIMA (2,0,2)(1,1,1)[12] | ARIMA (0,1,2)(0,1,1)[12] | |
|---|---|---|---|---|---|
| error model | 85.78 | 126.67 | 122.62 | 130.56 | |
| error predicció | 1977.984 | 218.4594 | 117.32 | 137.64 | 257.15 |